################################################################################
#################18.Robustness_alt_econ.R#######################################
################################################################################
#####Load plm package#####

library(plm)
library(pcse)
library(lmtest)
library(dplyr)
#Load Comparative Political Data Set

library(haven)
cpd <- read_dta("CPDS-1960-2016-Update-2018.dta") %>%
  select(year, country, openc, realgdpgr )

library(DataCombine)

#Lag variables and rename the UK#

cpd <-slide(data = cpd, Var = 'openc',TimeVar = 'year', GroupVar = 'country', slideBy = -1, NewVar = 'openc_lag')

cpd <-slide(data = cpd, Var = 'realgdpgr',TimeVar = 'year', GroupVar = 'country', slideBy = -1, NewVar = 'realgdpgr_lag')

cpd$country <- ifelse(cpd$country == 'United Kingdom', 'Great Britain', cpd$country)

#Rename country and year variables#

cpd <- cpd %>%
  rename(country.name = country,
         electionyear = year)

pols <- readRDS('pols_bjps')

#Merge the cpd data into pols#

pols <- left_join(pols, cpd, by = c('country.name', 'electionyear'))

#Scale openness and gdp growth#

pols$scale_open <- c(scale(pols$openc_lag))

pols$scale_gdp <- c(scale(pols$realgdpgr_lag))

######Run Model with openness instead of unemployment rate#####

pols <- pdata.frame(pols, index = c('country.name', 'order'))

full1 <- plm(scale_econsdw ~ scale_dispgini + 
               scale_ggini_inc1 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_ggini_inc1 + 
               scale_dispgini:scale_socioeconomicsal +
               scale_ggini_inc1: scale_socioeconomicsal +
               scale_ggini_inc1:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               openc_lag,
             data = pols, 
             na.action = na.omit, 
             model = "within")

coeftest(full1, vcov = function(x) vcovHC(x))

full2 <- plm(scale_econsdw ~ scale_dispgini + 
               scale_ggini_inc2 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_ggini_inc2 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_ggini_inc2: scale_socioeconomicsal +
               scale_ggini_inc2:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               openc_lag,
             data = pols, 
             na.action = na.omit, 
             model = "within")

coeftest(full2, vcov = vcovHC)

full3 <- plm(scale_econsdw ~ scale_dispgini + 
               scale_ggini_inc3 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_ggini_inc3 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_ggini_inc3: scale_socioeconomicsal +
               scale_ggini_inc3:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               openc_lag,
             data = pols, 
             na.action = na.omit, 
             model = "within")

coeftest(full3, vcov = vcovHC)

full4 <- plm(scale_econsdw ~ scale_dispgini + 
               scale_ggini_inc4 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_ggini_inc4 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_ggini_inc4: scale_socioeconomicsal +
               scale_ggini_inc4:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               openc_lag,
             data = pols, 
             na.action = na.omit, 
             model = "within")

coeftest(full4, vcov = vcovHC)

full5 <- plm(scale_econsdw ~ scale_dispgini + 
               scale_ggini_inc5 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_ggini_inc5 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_ggini_inc5: scale_socioeconomicsal +
               scale_ggini_inc5:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               openc_lag,
             data = pols, 
             na.action = na.omit, 
             model = "within")

coeftest(full5, vcov = vcovHC)

##### Table 1: Model 4 Output#####

fit1c <- c(summary(full1)$r.squared, summary(full1)$fstatistic$statistic, summary(full1)$fstatistic$p.value)
fit2c <- c(summary(full2)$r.squared, summary(full2)$fstatistic$statistic, summary(full2)$fstatistic$p.value)
fit3c <- c(summary(full3)$r.squared, summary(full3)$fstatistic$statistic, summary(full3)$fstatistic$p.value)
fit4c <- c(summary(full4)$r.squared, summary(full4)$fstatistic$statistic, summary(full4)$fstatistic$p.value)
fit5c <- c(summary(full5)$r.squared, summary(full5)$fstatistic$statistic, summary(full5)$fstatistic$p.value)


fit.statsc <- rbind(fit1c, fit2c, fit3c, fit4c, fit5c)

mean(fit.statsc[,1]) #R-Squared
mean(fit.statsc[,2]) #Adjusted R-Squred
mean(fit.statsc[,3]) #F-Statistic
mean(fit.statsc[,4]) #p-value of F-Statistic

fu1 <- as.data.frame(summary(full1, vcov = vcovHC)$coefficients) 

fu2 <- as.data.frame(summary(full2, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate2 = Estimate, 
                Std.Error2 = `Std. Error`, 
                `t-value2` = `t-value`, 
                `Pr(>|t|)2` = `Pr(>|t|)`)
fu3 <- as.data.frame(summary(full3, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate3 = Estimate, 
                Std.Error3 = `Std. Error`, 
                `t-value3` = `t-value`, 
                `Pr(>|t|)3` = `Pr(>|t|)`)
fu4 <- as.data.frame(summary(full4, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate4 = Estimate, 
                Std.Error4 = `Std. Error`, 
                `t-value4` = `t-value`, 
                `Pr(>|t|)4` = `Pr(>|t|)`)
fu5 <- as.data.frame(summary(full5, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate5 = Estimate, 
                Std.Error5 = `Std. Error`, 
                `t-value5` = `t-value`, 
                `Pr(>|t|)5` = `Pr(>|t|)`)

dot.data.fu <- cbind(fu1, fu2, fu3, fu4, fu5)

dot.data.fu <- dot.data.fu[ , order(names(dot.data.fu))]

dot.table.fu <- dot.data.fu %>%
  transmute(pe = rowMeans(dplyr::select(dot.data.fu, Estimate, Estimate2, Estimate3, Estimate4, Estimate5)),
            se = rowMeans(dplyr::select(dot.data.fu, `Std. Error`, Std.Error2, Std.Error3, Std.Error4, Std.Error4, Std.Error5)),
            t.value = rowMeans(dplyr::select(dot.data.fu, `t-value`, `t-value2`, `t-value3`, `t-value4`, `t-value5`)), 
            p.value = rowMeans(dplyr::select(dot.data.fu, `Pr(>|t|)`, `Pr(>|t|)2`, `Pr(>|t|)3`, `Pr(>|t|)4`, `Pr(>|t|)5`))) %>%
  mutate(lwr = pe - se*1.96) %>%
  mutate(up = pe + se*1.96)

dot.table.fu$Predictor <- c("Disposable Gini", "Income Dif.", "Economic Salience", "GALTAN Polarization", "ENEP", "Openness", "Gini*Income Dif.", "Gini*Salience", "Income Dif.*Salience", "Gini*Income Dif.*Saience")

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, 
                                 levels = c( "Income Dif.", "Economic Salience", "Disposable Gini","Income Dif.*Salience", "Gini*Income Dif.", "Gini*Salience","Gini*Income Dif.*Saience", "GALTAN Polarization","ENEP", "Openness"))

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, levels=rev(levels(dot.table.fu$Predictor)))

library(ggplot2)

#dot.table.fu %>% ggplot(aes(x = Predictor, y = pe)) + 
#  geom_point() + 
#  geom_hline(aes(yintercept = 0), linetype = 'dashed') +
#  geom_linerange(aes(ymin = lwr, ymax = up))+
#  coord_flip() +
#  labs(title = "Point Estimates of Model Regressing Predictors #on Economic Polarization (n = 83)", y = "Point Estimates (95% #Confidence Interval)", x = "Predictor")+
#  theme_bw()

rownames(dot.table.fu) <- dot.table.fu$Predictor
dot.table.fu <- dot.table.fu[,-7]

print(xtable::xtable(dot.table.fu, digits = 3))

#####Robustness check with gdp growth##### 


full1 <- plm(scale_econsdw ~ scale_dispgini + 
               scale_ggini_inc1 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_ggini_inc1 + 
               scale_dispgini:scale_socioeconomicsal +
               scale_ggini_inc1: scale_socioeconomicsal +
               scale_ggini_inc1:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_gdp,
             data = pols, 
             na.action = na.omit, 
             model = "within")

coeftest(full1, vcov = function(x) vcovHC(x))

full2 <- plm(scale_econsdw ~ scale_dispgini + 
               scale_ggini_inc2 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_ggini_inc2 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_ggini_inc2: scale_socioeconomicsal +
               scale_ggini_inc2:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_gdp,
             data = pols, 
             na.action = na.omit, 
             model = "within")

coeftest(full2, vcov = vcovHC)

full3 <- plm(scale_econsdw ~ scale_dispgini + 
               scale_ggini_inc3 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_ggini_inc3 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_ggini_inc3: scale_socioeconomicsal +
               scale_ggini_inc3:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_gdp,
             data = pols, 
             na.action = na.omit, 
             model = "within")

coeftest(full3, vcov = vcovHC)

full4 <- plm(scale_econsdw ~ scale_dispgini + 
               scale_ggini_inc4 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_ggini_inc4 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_ggini_inc4: scale_socioeconomicsal +
               scale_ggini_inc4:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_gdp,
             data = pols, 
             na.action = na.omit, 
             model = "within")

coeftest(full4, vcov = vcovHC)

full5 <- plm(scale_econsdw ~ scale_dispgini + 
               scale_ggini_inc5 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_ggini_inc5 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_ggini_inc5: scale_socioeconomicsal +
               scale_ggini_inc5:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_gdp,
             data = pols, 
             na.action = na.omit, 
             model = "within")

coeftest(full5, vcov = vcovHC)

##### Table 1: Model 4 Output#####

fit1c <- c(summary(full1)$r.squared, summary(full1)$fstatistic$statistic, summary(full1)$fstatistic$p.value)
fit2c <- c(summary(full2)$r.squared, summary(full2)$fstatistic$statistic, summary(full2)$fstatistic$p.value)
fit3c <- c(summary(full3)$r.squared, summary(full3)$fstatistic$statistic, summary(full3)$fstatistic$p.value)
fit4c <- c(summary(full4)$r.squared, summary(full4)$fstatistic$statistic, summary(full4)$fstatistic$p.value)
fit5c <- c(summary(full5)$r.squared, summary(full5)$fstatistic$statistic, summary(full5)$fstatistic$p.value)


fit.statsc <- rbind(fit1c, fit2c, fit3c, fit4c, fit5c)

mean(fit.statsc[,1]) #R-Squared
mean(fit.statsc[,2]) #Adjusted R-Squred
mean(fit.statsc[,3]) #F-Statistic
mean(fit.statsc[,4]) #p-value of F-Statistic

fu1 <- as.data.frame(summary(full1, vcov = vcovHC)$coefficients) 

fu2 <- as.data.frame(summary(full2, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate2 = Estimate, 
                Std.Error2 = `Std. Error`, 
                `t-value2` = `t-value`, 
                `Pr(>|t|)2` = `Pr(>|t|)`)
fu3 <- as.data.frame(summary(full3, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate3 = Estimate, 
                Std.Error3 = `Std. Error`, 
                `t-value3` = `t-value`, 
                `Pr(>|t|)3` = `Pr(>|t|)`)
fu4 <- as.data.frame(summary(full4, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate4 = Estimate, 
                Std.Error4 = `Std. Error`, 
                `t-value4` = `t-value`, 
                `Pr(>|t|)4` = `Pr(>|t|)`)
fu5 <- as.data.frame(summary(full5, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate5 = Estimate, 
                Std.Error5 = `Std. Error`, 
                `t-value5` = `t-value`, 
                `Pr(>|t|)5` = `Pr(>|t|)`)

dot.data.fu <- cbind(fu1, fu2, fu3, fu4, fu5)

dot.data.fu <- dot.data.fu[ , order(names(dot.data.fu))]

dot.table.fu <- dot.data.fu %>%
  transmute(pe = rowMeans(dplyr::select(dot.data.fu, Estimate, Estimate2, Estimate3, Estimate4, Estimate5)),
            se = rowMeans(dplyr::select(dot.data.fu, `Std. Error`, Std.Error2, Std.Error3, Std.Error4, Std.Error4, Std.Error5)),
            t.value = rowMeans(dplyr::select(dot.data.fu, `t-value`, `t-value2`, `t-value3`, `t-value4`, `t-value5`)), 
            p.value = rowMeans(dplyr::select(dot.data.fu, `Pr(>|t|)`, `Pr(>|t|)2`, `Pr(>|t|)3`, `Pr(>|t|)4`, `Pr(>|t|)5`))) %>%
  mutate(lwr = pe - se*1.96) %>%
  mutate(up = pe + se*1.96)

dot.table.fu$Predictor <- c("Disposable Gini", "Income Dif.", "Economic Salience", "GALTAN Polarization", "ENEP", "GDP Growth", "Gini*Income Dif.", "Gini*Salience", "Income Dif.*Salience", "Gini*Income Dif.*Saience")

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, 
                                 levels = c( "Income Dif.", "Economic Salience", "Disposable Gini","Income Dif.*Salience", "Gini*Income Dif.", "Gini*Salience","Gini*Income Dif.*Saience", "GALTAN Polarization","ENEP", "GDP Growth"))

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, levels=rev(levels(dot.table.fu$Predictor)))

library(ggplot2)

#dot.table.fu %>% ggplot(aes(x = Predictor, y = pe)) + 
#  geom_point() + 
#  geom_hline(aes(yintercept = 0), linetype = 'dashed') +
#  geom_linerange(aes(ymin = lwr, ymax = up))+
#  coord_flip() +
#  labs(title = "Point Estimates of Model Regressing Predictors #on Economic Polarization (n = 83)", y = "Point Estimates (95% #Confidence Interval)", x = "Predictor")+
#  theme_bw()

rownames(dot.table.fu) <- dot.table.fu$Predictor
dot.table.fu <- dot.table.fu[,-7]

print(xtable::xtable(dot.table.fu, digits = 3))